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We study here the steady state attained in a granular gas of inelastic rough spheres that is subject 
to a spatially uniform random volume force. The stochastic force has the form of the so-called white 
noise and acts by adding impulse to the particle translational velocities. We work out an analytical 
solution of the corresponding velocity distribution function from a Sonine polynomial expansion 
that displays energy non-equipartition between the translational and rotational modes, translational 
and rotational kurtoses, and translational-rotational velocity correlations. By comparison with a 
numerical solution of the Boltzmann kinetic equation (by means of the Direct Simulation Monte 
Carlo method) we show that our analytical solution provides a good description that is quantitatively 
very accurate in certain ranges of inelasticity and roughness. We also find three important features 
that make the forced granular gas steady state very different from the homogeneous cooling state 
(attained by an unforced granular gas). First, the marginal velocity distributions are always close 
to a Maxwellian. Second, there is a continuous transition to the purely smooth limit (where the 
effects of particle rotations are ignored). And third, the angular translational-rotational velocity 
correlations show a preference for a quasiperpendicular mutual orientation (which is called “lifted- 
tennis-ball” behavior). 


I. INTRODUCTION 

Systems composed of a large number of particles with a characteristic length bigger than about 1 pm 
(granular systems) are ubiquitous in nature^ and applications^ Therefore, the study of the transport 
processes in this kind of systems (generically called granular matter ) has been of interest and subject of 
systematic study for a long time (see, for instance, the work by Reynolds^). The most generic property 
of granular matter is the non-conservation of kinetic energy after contact between particles, and hence the 
term inelastic particles*^ 

Attending to the nature of the mechanical interaction among grains, granular matter may be classified 
into two groups: wet granular matter (grains tend to stick together after contact and they are usually called 
cohesive powders , since they usually predominate in the smaller range of grain sizes^) and dry granular 
matter (grains do no not show this sticky behavior). Of course, both types of granular systems are equally 
relevant in industry and technologyjii^&i Moreover, the finding of a meaningful theoretical description also 
poses fundamental challenges for researchers of a number of fields^ Since the nature of particle contacts 
essentially determines the granular dynamics, theoretical developments frequently study granular systems 
as composed by only dry or only wet particles (or, at its most elaborate, simple combinations of both). 

Let us focus on dry granular matter. Here, we can find persistent contacts as well, but they are always 
due to a high fraction of the volume occupied by grains. When this happens, available space limitations 
result in lasting contacts. As a consequence, mechanical force transmission throughout the system occurs*^ 
i.e., one particle may be in contact with not just one another particle but a set of them, thus propagating 
mechanical information to large parts of the system^ Correlations in this case are obviously strong, much in 
the same way as in atomic solid matter^ On the contrary, in less dense systems, particles are less correlated 
and contacts among three or more grains tend to be statistically rare events; i.e., contacts typically become 
only binary^ At the same time, in the low density regime, those contacts between pairs of dry grains tend 
to be very short in time (as compared with the characteristic mean free time), and hence the term inelastic 
collision d 1 Binary and instantaneous collisions are also characteristic features of low density molecular gases, 
for which the so-called molecular chaos assumption ( Stosszahlansatz ) is commonly usedji 2,13 i.e., colliding 
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particle velocities are supposed to be uncorrelated. In those circumstances, the Boltzmann equation or, at 
low but finite density, the Enskog equation (properly modified in order to take into account inelasticity in 
the collision s 11 ! 13 ) applies. While inelastic collisions may originate velocity correlations,— these effects tend 
to be more important at higher densities^ and thus they can usually be ignored for low densities/ ' 11 ' 16 Since 
the theoretical approach to both the high and low density limits for granular matter is radically different/ 
(actually, just in the same way as it is for conventional matter—), it is convenient to focus on just one density 
regime. In this case we are interested in the low density limit, the so-called granular gases.— 

Given that most granular particles do not have a regular shape and differ in size and mass/^ granular 
collisions (even if short in time) can be very complex events ! 20 ' 21 It is thus important that granular colli¬ 
sion models can accordingly capture as many features as possible, without otherwise compromising their 
tractability. For instance, the coefficient of normal restitution (that characterizes translational kinetic en¬ 
ergy loss upon collisions) depends in general on the impact velocit y 11 ' 21 (a viscoelastic collision model may 
be used to take this into account—). Nevertheless, it has been found that an important set of the main 
features at a macroscopic level (i.e., average field behavior) can be described with models that ignore the 
effects of impact velocity on the restitution properties.- In an effort to simplify the collision process, one 
can sequentially describe it as (i) an impact between both colliding particles, after which they suffer a me¬ 
chanical impulse that produces a linear momentum transfer; (ii) a friction process, tangent to the contact 
interface, generating an angular momentum transfer; and (iii) a sliding process, that describes the relative 
translational displacement of the colliding particles in the direction of the contact interface. All these three 
features were captured in previous research, that includes experimental work (see for instance Ref. [23h , by 
using a collision model based on three collisional constant parameters: a coefficient of normal restitution, a 
coefficient of tangential restitution, and a friction coefficient. 

In spite of the intrinsic polydispersity of granular systems in nature/ a variety of generic granular dynamics 
phenomena can be described by means of kinetic theory for a monodisperse granular gas. In a granular 
gas of identical particles, the relevance of collision models based on constant coefficients of restitution 
is exemplified by the correct description that the simple inelastic smooth hard-sphere model gives of the 
clustering instability of the so-called homogeneous cooling state (HCS), i.e., an unforced state with uniform 
hydrodynamic quantities and a time decaying granular temperature/ 4 ~— Furthermore, simplifications of the 
collision model alleviate the involved mathematical task of integrating the Boltzmann equation, especially for 
states more complex than the HCSi 27 ' 28 In turn, solutions of the Boltzmann equation have allowed for finding 
new fundamental properties of granular gas dynamics by means of consistent derivations of the corresponding 
hydrodynamic theories— Another important example of this is the inelastic rough hard-sphere model, that is 
able to detect correlations between the angular and translational velocity components in granular gases . 29 ' 30 
Therefore, it is fruitful to compromise between simple collision models and more sophisticated theoretical 
developments at a kinetic theory level, in order to discover new fundamental transport properties of granular 
gases. 

In its basic implementation, the rough hard-sphere model neglects the effects of eventual sliding in collisions 
but is capable of showing a number of relevant results for the dynamics of granular gases. For instance, 
it has shown that energy non-equipartition between translational and rotational modes occurs/* — and 
that the distribution function of particle velocities can exhibit strong non-Maxwellian effects—— It has 
also been determined that strong inelasticity or angular-translational correlations do not necessarily lead to 
difficulties for reaching a normal (or hydrodynamic) state. 36 This observation has important consequences for 
the development of hydrodynamic theories for granular gases/’ 3 '. However, most of these works addressed 
only the HCS. Useful and fundamental as this state may be as a reference for transport theories of granular 
gases//2Z it is hard to observe its spontaneous appearance in nature and also difficult to reproduce in 
laboratory experiments. In fact, granular gases are most frequently observed under a forcing that keeps the 
system steadily fluidized and at low density. 

In this paper, we study the homogeneous steady state of a granular gas of rough spheres excited by 
a stochastic volume force that is modeled by means of the so-called white noise^r— as usually named in 
statistical physics. We extend previous theoretical results^ by taking into account the effects of orientational 
correlations in particle velocities and, in addition, by comparing the theoretical predictions with computer 
simulations. As we will see, even though the velocity kurtoses and correlations are still not negligible, 
the heated granular gas has a more Maxwellian-like distribution function for all values of the roughness 
parameter than the HCSi 30 ' 36 This justifies the practical validity of our theoretical perturbation approach 
(truncated Sonine expansion) and its excellent agreement with simulation data. In contrast, the HCS is 
characterized by large deviations from the Maxwellian distribution for small roughness (especially in what 
concerns the angular velocity), in which case the perturbation method is valid only semi-quantitatively— 
Moreover, no singularity exists in the smooth-sphere limit for the heated gas, again in contrast to what 
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happens in the HCS. 

Another important and distinctive feature that we will find is that the forced granular gas shows an average 
preference for quasiperpendicular relative orientation of translational and rotational particle velocities, in 
contrast to the behavior found in the unforced granular gas , 29 ' 30 ' 36 ' 44 where quasiparallcl translational and 
rotational particle velocities are also present. This result is of interest for a number of experimental works, 
where the most common situation is that of an excited granular gas (see, for instance, Ref. HU). 

The structure of this work is as follows. In Sec. El we write the corresponding kinetic (Boltzmann) 
equation. By defining adequate time and velocity units, we also derive the dimensionless counterparts of 
the Boltzmann equation and its derived moment equations, which are the ones we will work with. In Sec. 
nm a theoretical solution to the moment equations, by means of a Sonine polynomial expansion around the 
Maxwellian distribution function, is derived by keeping terms up to order two (fourth-degree polynomials). 
In this way, the state of the system is solved at all times , including the final steady state. In particular, the 
theoretical parameters characterizing kurtoses, correlations, and average relative orientation of translational 
and rotational particle velocities are analyzed. As said before, we will find that excited rough grains do not 
show the relevance of quasiparallcl relative orientation behavior present in the HCS . 29 ’ 36 Next, in Sec. IIVI 
we compare the theoretical results with an “exact” numerical solution of the Boltzmann equation obtained 
by means of the Direct Simulation Monte Carlo (DSMC) method,— a very good agreement being observed. 
Finally, the results are briefly discussed in Sec. |V] 


II. BOLTZMANN EQUATION AND MOMENT HIERARCHY 
A. Collision rules 


Let us consider a homogeneous dilute granular gas of identical hard spheres with mass to, diameter a, and 
moment of inertia /, subject to a stochastic volume force F wn (also called thermostat) with the properties 
of a Gaussian white noise— ~ 42 i 47 This kind of forcing can model, for example, the energy input to grains 
immersed in a gas in turbulent flow— 

We denote the translational and angular (or rotational) particle velocities with v and uj respectively. 
Binary collisions are characterized here by two material parameters, namely the coefficient of normal resti¬ 
tution a and the coefficient of tangential restitution /3, which determine the shrinking of the normal and 
tangential component, respectively, of the relative velocity of the two surface points at contact. They are 
defined by (see, for instance, Refs. [37l lilj. and [Hdl 

a ■ u' = — a(a ■ u), a x u' = — /3(<r x u). (2-1) 


Here, the primes denote post-collisional values, <r is the unit collision vector joining the centers of the two 
colliding spheres (and pointing from the center of particle 1 to the center of particle 2), and 


U = Vi - v 2 - —a x (uq + uj 2 ) 


( 2 . 2 ) 


is the relative velocity of the points of the spheres which are in contact during a binary encounter. 

The coefficient of normal restitution a takes values between 0 (completely inelastic collision) and 1 (com¬ 
pletely elastic collision), while the coefficient of tangential restitution takes values between —1 (completely 
smooth collision, implying that rotational velocities are unchanged) and 1 (completely rough collision) d 2 ' 49 
Equation supplements the linear and angular momentum conservation laws to yield the collision 

rules, 11 ' 35 ' 37 ' 49 


mV 'l,2 = mvi-2 T Q: lu 1,2 = 1^1,2 - -CT X Q, 

where the impulse exerted by particle 1 on particle 2 is given by 

Q = ma(<r ■ u)<x — m/3<r x (<r x u). 

Here, we have introduced the abbreviations 


1 T cr 'r _ 1 T /? Si 

1 P — 


4 1 


t2 


(2.3) 


(2.4) 


a = 


2 


2 «+ 1’ 


K — 


(2.5) 
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While collisions (except if a = 1 and /3 = ±1) dissipate kinetic energy (translational plus rotational), 
the homogeneous stochastic force F wn injects translational kinetic energy to grains. It has properties of a 
Gaussian white noise, i.e., 


(F 7 n W) =0, (2.6a) 

(FrWFr(O) =\m 2 xlS Z] S(t - t 1 ), (2.6b) 

where indexes i, j refer to particles, I is the 3x3 unit matrix, and Xo measures the characteristic strength 
of the stochastic force. 


B. Boltzmann equation 


In homogeneous states, the Boltzmann equation corresponding to the stochastic external force F wn 
become s 41 ’ 42 


df(v,u;t) 

dt 


Xo 



2 

/(v,u>;t) = J[v,w|/(f)], 


(2.7) 


where /(v, u>;i) is the velocity distribution function and J[v,u;|/] is the collision integral in the (inelastic) 
Boltzmann equation for rough spheres, which accounts for the collision rules (12.311 and (|2.4I) . 

Given a certain function H(v,u;), we define its average as 

(A(t)) = i J dw J rfwA(v,w)/(v,«;t), n = J dv J doj /(v,u>;f), (2.8) 

where n is the number density. By Galilean invariance, we choose (v) = 0. Moreover, we assume isotropic 
states, so that (<*>) = 0. 37 Thus, the basic quantities are the translational (T t ), rotational (T r ), and total 
(T) granular temperatures: 



Tr 



_T t + T r _ 1 + 6 

2 4 2 


(2.9) 


where the temperature ratio 9 = T r /T t is an important quantity in this system because, as we will see, its 
steady-state value is independent of the level of forcing (xg)- 

The evolution equations for the granular temperatures T t and T r are obtained by just multiplying both 
sides of Eq. m by the particle translational and rotational kinetic energies, respectively, and integrating 
over all velocity values. This gives 


dtTt — m \q — — £tT t , d t T r — —t; r T r , d t T — 


mXo 


= -C T 


( 2 . 10 ) 


Here, the translational (£ t ) and rotational (^ r ) production rates and the cooling rate (£) are defined by the 
collision integral s 35 ! 37 


it = ~ 

c = 


JdvJ dui v 2 J[v,u:\f], Cr = -^r J dv J duj w 2 J[v, u>\f], 


3 nT t 

£tTt + £, r T r it + C, 


2 T 


i + e 


(2.11a) 

(2.11b) 


C. Reduced variables and velocity moments 


For convenience, let us rewrite the Boltzmann equation m in dimensionless form. For this, we first 
define the reduced velocity distribution function 


</>(c,w ;r) 


1 ^ 4 T t {t)T r (t) \ 3/2 


ml 


/(v,w;t) 


( 2 . 12 ) 
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where the reduced particle velocities are 

V bJ 

C (t) = =. W (t) = =. 

y/2T t (t)/m V2 T r (t)/I 

In addition, time is measured by the parameter r defined by 

r(t) = ( dt'v{t'), v{t) = 2na 2 \JTrT t (t) / m, 

Jo 


(2.13) 


(2.14) 


where v(t) is the (time-dependent) collision frequency. Thus, r(f) measures the accumulated number of 
collisions per particle up to time t. With this time scale, we can measure relaxation times to the steady 
state that are physically meaningful, since they coincide with the aging times to hydrodynamics^ 

The moments of the reduced distribution function c, w; r) are 


Mpq\r) = (c p w q (c ■ w) r ) 



dw c p w q (c ■ w) r </>(c, w; r), 


(2.15) 


(r) . (0) 

with p,q,r = even, by symmetry. Note that Mp q is a moment of degree p + q + 2r and Mg J = 1, 
M^o = Mq °2 = • The dimensionless measure of the noise intensity is defined by 


= 3 Xo = _^Xo_ 

2 v(t)T t (t)/m 4:^/nna 2 [T t (t)/m\ 3 / 2 


(2.16) 


Notice that, although Xo is a constant, its dimensionless counterpart 7 is time-dependent due to the scaling 
with the temperature T t . Of course, once a steady state is reached, 7 becomes constant too. Note also 
that 7“ 2 / 3 can be seen as the translational temperature Tt in units of a reference temperature T re f = 
m(3xo/4:y/Trna 2 ) 2 / 3 defined from the white noise parameter Xo- 
Let us finally define the reduced collisional moments 


where 



dw c p w q (c • 


w) r J[c,w|()>(T)], 


J[c,-w\(j){T)\ = 


AT t (t)T r (t) 


(t) |_ ml 


3/2 


■J[v,w|/(t)] 


(2.17) 


(2.18) 


is the dimensionless form of the collision integral. From the definitions (12.llal) and (12.171) . it is easy to see 

that | u IJ-w ; Cr = 'i u l 4)2 ■ The evolution equations for the temperature ratio 6{t) and the reduced noise 

strength y(r) can be obtained from Eq. (I2.10[) as 


d hi 6(t) 
dr 


2 

3 


a4o( t ) 


M 02 ( t ) ~ 7 (t) 


9 In 7 (r) 


M 2 o( r ) 


(2.19) 


Note that, since v and u) are scaled each with a different temperature in Eq. (12.131) . the collision rules (12.31) 
and (G3D, when expressed in terms of c and w, depend parametrically on 9 35 This dependence is obviously 
transferred to the collisional moments (12.171) . 

Taking into account the previous dimensionless quantities, and making use again of the temperature 
evolution equations (|2.10D . one can obtain the dimensionless version of the Boltzmann equation (12.71) as 


9(/>(c,w;r) p4o( t ) "TOO d 
dr 3 dc 


M>(c, w; t)] + 


M02V) 9 

3 dw 


7(t) 

6 



2 

0(c,w;t) = J[c,w\(j){T)\. 


W(c,w;r)] 


( 2 . 20 ) 


Multiplying both sides by c p w q (c ■ w) r and integrating over c and w, we get the moment hierarchy 
91nMp^(r) p + r r — n "- 1 -- — 


dr 


/4o ( T ) - 7 (t) + - 4 —^02 ; ( r ) - 


q + r ( 0 )/ ^ V P q (r) 




7 (r) P(P+ 1 + 2 r)M^ q (r) + r(r - l)M^ g+ 2 2 } (r) 
Mpq ( T ) 


6 


( 2 . 21 ) 
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This equation must be supplemented with Eq. (12.191) . As will be seen later, the steady-state distribution 
function is independent of the precise value of the stochastic force intensity Xo or > equivalently, of the initial 
value 7 ( 0 ), although of course the transient states are dependent on it. 


D. Exact Sonine expansion 


In principle, the reduced distribution function <j>(c, w) would be a function of six velocity components. 
On the other hand, since we restrict ourselves to isotropic states, <f> is actually a function of three scalar 
quantities: the norms c 2 and w) 2 , and the squared dot product (c, w) 2 ^ Thus, it makes sense to represent 
(j){ c,w) as an expansion in a complete set {'I'^ (c, w)} of orthogonal polynomials, 


OO OO OO 


0(c,w;r) = <j> M (c,w) = tt 3 e 

j =0 *;=0 e=o 


—6„—C —w 


where 


*$( c-w 


) = L^ t+ 2} (c 2 )L^ e+ * ) (w 2 ) (c 2 w 2 Y P 2f .{cosd), cos'd = 


cw 


( 2 . 22 ) 


(2.23) 


where ZA 2 (x) and P 2 e(x) are the Laguerre and Legendre polynomials, respectively^ Defining the scalar 


product of two arbitrary isotropic (real) functions $i(c, w) and <I> 2 (c, w) as 

($i|$ 2 ) = J dc j dw 4> m {c, u;)$i(c, w)<h 2 (c, w), 


one can obtain the orthogonality relation 

= Mjkhi'8k,k>S t s, = 

As a consequence, the coefficients in the expansion (12.221) are 


r(j + 2 £+|)r(fc + 2 £ + |) 

( 7 r/ 4 )( 4 £ + l)j!fc! 


(2.24) 


(2.25) 




= 


( l )jl+fcl+ , l(7r/ 4 )j!fc!(4£ + 1)(4£ l)W%££$ Xki+li 


££>£ P'-ti-JiW (ii +2^+1) h\(k - h)\r (k! + 2 £+l) 2 ‘i£ 1 l( 2 £- 24)!’ 


(2.26) 


where in the second step use has been made of the explicit expressions of the Laguerre and Legendre 
polynomials^ We observe that a'^ is a linear combination of the moments Mp T q with p 7 q 7 r = even and 

degree 21 < p + q + 2r < 2 (j + k + 21). By normalization, Oqq* = 1, a\ = a ^ = 0. Therefore, the first 
nontrivial coefficients are those associated with moments of degree four, namely 

4 


4o = ^|(c 4 ) - 1, a ( 0 °2 = 44 w ^ ~ 1} 


(2.27a) 


jo) - — / c 2 w 2y > - 1 a (1) - A 

a ll ~ g\ C W / A > a 00 ~ 25 


((c • w) 2 } — —(c 2 w 2 ) 


(2.27b) 


These are the fourth-degree cumulants , which measure the basic deviations of the velocity distribution 
function w) from the (two-temperaure) Maxwellian (/>m(c,w). While and measure the kurtosis 
of the (marginal) translational and rotational velocity distribution functions, respectively, the coefficient 
a^i characterizes the scalar translational-rotational correlations, i.e., (c 2 w 2 ) J (c 2 )(w 2 ). Finally, a ^ 
informs us about the possible existence of orientational correlations, i.e., ( c 2 w 2 cos 2 d) J (c 2 w 2 )( cos 2 d) and 
(cos 2 d) J To disentangle the latter two effects, let us define the quantities 




((c-w) 2 ) 


(c 2 w 2 )( cos 2 d) 


- 1 


6 = W ( ( cos2 #> - 7 


(2.28) 
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The parameters h, b , a 


(o) 
n i 


and a 


(i) 

oo 


are related by 


1 25a$ - 96 (l + 

16 l 1 + To b ) ( x + a u) 


(2.29) 


In previous works ) 29 ’ 30 ' 44 the angle •& between translational and rotational particle velocities has been used 
for measuring orientational correlations in a granular gas of rough spheres, by tracking the quantity (cos 2 d) 
or, equivalently, b. The use of b has the advantage over a qq* of being more intuitive (since it directly refers 
to the angle between translational and rotational particle velocities). However, it has the disadvantage that, 
since cos'd = c • w/cw is not a polynomial, the parameter b involves an infinite number of coefficients in the 
expansion (12.221) . More specifically, it is easy to find 


OO OO 

6 =EE a $- ( 2 - 3 °) 

j —0 k =0 

Analogously, the quantity h also involves an infinite number of coefficients. For this reason, we prefer 
instead to use the single coefficient for measuring the orientational correlations in a granular gas of 
rough spheres^ 

Taking into account the general moment hierarchy (12.21 II . the evolution equations for the fourth-degree 
cumulants are 


d In 


1 + <4o 


0) 


dr 


4 (0) 
= 3^° 


(t) _ iy( r) - AwiElM 
( ) 3 " ( ) 15 i + 4°oE) 


(2.31a) 


d In 


1 + a 02 ( T ) 


J _ 4 ..( 0 ) 4 M04( t ) 


dr 


= 3^02(r)-- (Q) 

d 4t) 1 + % 2 \ T ) 


(2.31b) 


d In 


1 + a H ( T ) 


dr 


2 (o), \ , 2 (o), , 2 4/42 } (t) - §7(t) 

= 3^ 20 (T) + 3^ 02 (T) - 3 7(T) “ 9 i + a (0) (T) * 


(2.31c) 


d In 


1 + a u ( T ) + | a oo ) ( r ) 


dr 


2 (0). X , 2 (0), , 2 , ^ 4 MooE) - h( T ) 

= 3^20 M + 3M02 M - p ( r ) - 3 1 + a (0) (r) + |a (i) (r) - 


In the steady state (d T 0), Eqs. (12.101) and (12.311) imply the conditions 


_ (o) (o) _ n 

7 — d20 ) do2 — 0; 


(2.31d) 


(2.32a) 


I u (0) _ 2 (0) _ 2 (2) _ (0) (o)_ Q 

gd40 ~~ 3P22 — 00 — d20 1 P04 — u - 


(2.32b) 


Before closing this section, and for further use, let us define the marginal distributions </> c (c), (f> w (w ), 
cw(c 2 w 2 ), </>c w((c • w) 2 ), and ^(cos 2 d) from the joint velocity distribution function <(>(c, w) as 


4> c {c) = 47tc 2 J dw</>(c,w), <j> w (w) = 4nw 2 J dc</>(c,w), 


(2.33a) 


4>cw(c 2 w 2 ) = J dc 'j dw' i(c' 2 w' 2 - c 2 w 2 )()(c', w'), 


(2.33b) 




























0 c w((c • w) 2 ) = J dc' J dW S((c' ■ w ') 2 - (c • w) 2 )</>(c', w') 


(2.33c) 


tj>s( cos 2 i?) = J dc' J dW <5(cos 2 d' — cos 2 $)</>(c', w'). 


(2.33d) 


The less conventional marginal distribution functions (12.33bl) (12.33dll are directly related to the scalar and 
orientational correlations between the translational and rotational particle velocities. More specifically, 

p OO pO O pi 

(c 2 w 2 ) = / dxxcj) cw (x), ((c • w) 2 ) = / dxx(j> c .- w (x), (cos 2 d) = / dxx(j>#(x). (2-34) 

Jo Jo Jo 


III. APPROXIMATE SOLUTION FROM A TRUNCATED SONINE EXPANSION 

All of the equations in Sec. [Til are formally exact. On the other hand, the collisional moments 
1*1*02 , a4cm an< T l*oo > defined by the general expression 1)2.1711 . are functionals of <j>(c,w), and thus 

they depend on the whole set of expansion coefficient {a^}. They also depend on the temperature ratio 9 
through the collision rules. Therefore, Eqs. (12.191) and (12.311) do not make a closed set of equations. 

Nevertheless, the expansion ( 12 . 221 ) is especially useful for states with a distribution function close to the 
Maxwellian 4>m(c, w ), in which case an approximate closure can be applied. In a previous work^ we used 
the expansion ( 12 . 221 ) for describing the HCS distribution function. We found that the HCS is very far from 
the Maxwellian for a certain range of values of the roughness parameter 0, thus limiting the practical validity 
of this kind of theoretical description for the unforced gas. However, we can expect that the white noise 
thermostat keeps the distribution function closer to the Maxwellian, in analogy with what happens in the 
smooth-sphere case.— If so, a better accuracy of the expansion (|2.22l) can be expected in the present case. 
We will confirm this property in Sec. IIVI 


A. Maxwellian approximation 


The simplest approximation consists in neglecting all the cumulants , j + k + 2i> 2, in Eq. (12.221) . i.e., 
<j){ c, w) ~ <Pm{c,w). This Maxwellian approximation is of course unable to account for the evolution and 
steady-state values of the cumulants. However, it can capture the main features of the temperature ratio 
9{t) and the reduced noise intensity y(r). When </>(c, w) ~ is inserted into Eq. (12.171) . one obtains 

/4o!m = 1 — a 2 + (j + J ^2 P + K (-*- ~ P) ~ + P)\ > (3.1a) 

/$,m = [2 + « _1 (1 -P)~ 0 _1 (1 + P)\ ■ (3-lb) 

Using these expressions in (12.191) . one gets a closed set of two coupled nonlinear differential equations that 
can be numerically solved for arbitrary initial values 9(0) and q(0). Regardless of the values of 0(0) and 
7 (0), common steady-state values are reached after a few collisions per particle. They are explicitly obtained 
from Eq. (12.32al) as 


9 m = 


1+P 

2 + k -1 (1 — f3) ’ 


7m = 1 - a 2 + 


2(1 ~~ /3 2 ) 

2 + k ~ 1 (1 — ( 3 )' 


(3.2) 


It is interesting to note that the Maxwellian prediction for the steady-state value of the temperature ratio 
is independent of the coefficient of normal restitution a. Moreover, 9m is always smaller than unity (i.e., 
T r < T t ), except in the limit of completely rough spheres (0 = 1), where 9m = 1. In the case of the reduced 
noise intensity, we can observe that 7 m monotonically decreases as a increases at fixed 0, as a consequence 
of a progressively smaller cooling effect. Because of the same reason, 7 m presents a non-mo notonic behavior 
with respect to 0 at fixed a, reaching a maximum value 7 m, max = 1 — a 2 + 4k(1 + 2k — 2y/ k + k 2 ) at 
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0 = 1 — 2(y /k + k 2 — k). However, as we will see in Sec. hyi e actually depends on a but much more weakly 
than it depends on 0. Also, 6 can reach values slightly larger than unity if 0 is very close to 1 and a is small 
enough. In any case, as shown later, the Maxwellian expressions m are indeed rather accurate. 

The marginal distribution functions defined by (|2.33l) adopt simple forms in the Maxwellian approximation, 


<^c,m(c) 






(3.3a) 


0cw,m(c 2 w 2 ) = —— Ao(2cw;), </>c w,m((c • w) 2 ) = -K x { 2|c • w|), <j># ,m(cos 2 •&) = 1 , (3.3b) 

7r 7r 2 |cost7| 

where I\ n (x ) is a modified Bessel function of the second kindj^i Upon deriving (I3.3bl) use has been made 
of the mathematical properties S(c' 2 w' 2 — x) = S(w' — y/x/c')/2c ,2 tu\ S((c' • w') 2 — x) = <5(| cos i9| — 
\/x/cw)/2c 2 w 2 \ cos$|, and (5(cos 2 d — x) = <5(| cos$| — yfx)/2\ cos$|, respectively. Apart from (c 4 ) = (w; 4 ) = 
it is straightforward to check from Eqs. (12.341) and (13.3bl) that (c 2 w 2 ) = |, ((c-w) 2 ) = |, and (cos 2 $) = | 
in the Maxwellian approximation, as expected. 


B. Truncated Sonine approximation 


As a much more elaborate approximation that captures the most important non-Maxwellian features of 
the velocity distribution function, we truncate the exact expansion (12.221) after j + k + 21 = 2, i.e»2S, 

-l+^cMo (c, w) + ( c , w ) + «Mi ( c > w ) + a oo *oo ( c > w ) 

(PM{C,IVJ 


— 1 ~h CL , 


(0) 15 - 20c 2 + 4c 4 ( 0 ) 15 - 20w 2 + 4w 4 


20 


J '02 


( 0 ) 

*11 


(3 - 2c 2 ) (3 - 2 w 2 ) 


(i) 3(c • w) 2 — c 2 w 2 


+ CL , 


00 


(3.4) 


This approximation differs from that considered in Ref. |35| by the addition of the term headed by the 
coefficient . The associated marginal distributions are 


4> c {c) (o) 15 - 20c 2 + 4c 4 

, , \ — 1 + a 20 o > 

<Pc,M\ c ) 8 


4>w (w) (o) 15 - 20u> 2 + 4 w 4 

- - - = 1 + dg 2 -5-, 

<i Vw,MV W ) 8 


(3.5a) 


<j) cw (c 2 w 2 ) =1 + a^o + ay? 


4>cw,m(c 2 W 2 ) 


15 + 4c 2 u> 2 — 16cw 


Ki(2cw) 


A'o(2 c 


+ 


,(°) 

ni 


o r, K 1 (2 cw) 

9 + 4c 2 u> 2 — 12cu>' 


K 0 (2cw) 


(3.5b) 


</>c- w ((c ■ w) 2 ) 

</>c-w,m((c • w) 2 ) 


,(o) , _(o) 


=1 + 


^20 


^02 


3 + 4(c • w) 2 — 12|c • w 


A' 0 (2|c • w|) 


+ 


z (0) 

ni 


1 + 4(c • w) 2 — 8|c ■ w 


A'i(2|c • w|) 
A' 0 (2|c • w|) 


"f- 1 — 2(c • w) 2 + |c ■ 


A'i(2|c • w|)_ 

I A' 0 (2|c ■ w|) 

Ad(2|c • w|) 


(3.5c) 


COS 2 $) 

</VM (cos 2 l}) 


9a, 


(i) 


= 1 + —(3 cos 2 1 ? — 1) . 


(3.5d) 


This yields (c 4 ) = x(l+° 2 o), ( w ‘ 4 ) = T( 1 + a 02 ) ) I (c 2 ™ 2 ) = f( 1 + a u)> (( c ’ w ) 2 ) = |0 + «n + § a oo)> and 
(cos 2 tf) = 1(1 + j^a gg^), as expected. The latter equality implies b = a^, in consistency with the neglect 
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of a^/j with j + k + 2£> 3 in Eq. (12.301) . If, additionally, terms nonlinear in the cumulants are neglected in 
Eq. (12.291) . we have 


h ~ b ~ ago (3-6) 

in the Sonine approximation (13.41) . Thus, if < 0 (as will be seen to be the case for most values of a and 
/? in the case of uniform spheres), we can expect that (cos 2 $) < i. This implies a tendency of the vectors 
c and w to favor quasinormal mutual orientations (“lifted-tennis-ball” effect). 

If Eq. (13.41) is inserted into Eq. (12.171) and terms nonlinear in a^ k ! are neglected, one obtains explicit (yet 

not exact) expressions of /Jpq with p + q + 2r = 2 and 4. Those expressions were first displayed in the 
Supplemental material of Ref. [36] but, for the sake of completion, they can be found here in the Appendix. 
Once Eqs. (IA.1I) (IA. 6 I) are introduced into Eqs. (12.191) and (12.311) . the resulting closed set of six equations 
can be numerically solved for arbitrary initial values to get the time-dependent functions 6{r), y(r), a^(r), 

a 02 ( T )i anV), and ffloo^C 7 ’)'^ We include in Sec. [IV] graphs of those functions and compare them with 
the “exact” numerical solution of the Boltzmann equation obtained by means of the DSMC method^ 

As for the steady state, by applying the state-state conditions (12.321) and the Sonine expressions (IA. 1 I) 
(IA. 6 I) . we are able to obtain fully analytical explicit expressions for the temperature ratio 9 , the dimensionless 

noise intensity 7 , and the four cumulants a^o*, a 02 > a ii > and a cw- The method is quite simple^ First, the 
relations in Eq. (12.32bl) are used to express the four cumulants in terms of the yet unknown quantity 9. 
Inserting those expressions into pjyj = 0, one gets a closed quartic equation for 9, whose physical solution 
is identified as the one closer to 9m [see Eq. mi Finally, 7 is obtained from 7 = /jl^q ■ 

Figure Q] presents density plots of the steady-state values for uniform spheres (k = §) of 9, 7 , , ajy], 

, and cigg*. In Fig. [lja) we may point out that iso-0 curves are almost perpendicular to the /? axis, 
indicating that 9 depends almost entirely on the coefficient of tangential restitution and hardly on a, as 
anticipated by the Maxwellian approximation 9m hr Eq. (13.21) . Also, 9 < 1, except, and this is not visible in 
Fig. [2(a), if 0 > 0.994 336 and 0 < a < aq(/3), where aq(0.994 336) = 0 and aq(l) = l/\/2. The maximum 
value 0 max = 1.009 85 takes place at a = 0 and /3 = 1, i.e., in the completely inelastic and rough limit. 
The behavior of 7 observed in Fig. [2b) is in good agreement with the one expected from the Maxwellian 
approximation 7 m in Eq. & The latter predicts, at fixed a, a maximum value 7 m, max — 1-4853 — a 2 
at fj ~ 0.3033. In the Sonine approximation, we find that the location and value of the maximum slightly 
changes from /? ~ 0.301 and 7 max + a 2 — 1.512 at a = 0 to j3 ~ 0.304 and 7 max + a 2 — 1.486 at a = 1. 

Since 9 and 7 are well represented by the Maxwellian approximation (13.21) . the most interesting panels in 
Fig. [2 are those related to the cumulants. The first feature to be noted is that all of them have a relatively 
small magnitude |o^?| < 0.1, implying that a good performance of the Sonine approximation (13. 41) . and 
hence the validity of (13.61) . can be expected. This expectation will be confirmed in Sec. El We can also see 
that, at given a, the maximum value of |a^| is typically reached in the region of intermediate roughness 
(/? ~ 0), which differs from the behavior detected for the HCS, where the maximum values (which are also 
much larger than here) occur in the region close to the smooth limit .4 6 

Figure [2(c) shows that the translational kurtosis a^f takes negative values in the top and bottom right 
corners of the map. In contrast, the region where the rotational kurtosis is negative reduces to a small 
top-right corner (a,/3 < 1) [see Fig. [2d)]. The analogous top-right corners where < 0 and > 0, 
respectively, are even smaller in the cases of the correlation cumulants and, especially, [ see Figs. 
□(e) and [Q(f)]. 

Therefore, the typical non-Maxwellian features of the joint velocity distribution function cj)( c, w) that can 
be extracted from Figs. mcH 2 f) are: (i) platykurtic translational (d^l > 0 => (c 4 ) > ^) and rotational 
( a 02 > 0 =£* (u; 4 ) > velocity distributions; (ii) particles with larger (smaller) translational velocities 
tend to have larger (smaller) rotational velocities ( 07 ° > 0 =£* ( c 2 w 2 ) > (c 2 )(w 2 ) = |); and (iii) quasinormal 
orientations between the vectors c and w tend to be favored against quasiparallel orientations ( 0 ^ < 0 => 
(cos 2 iD) < i). Exceptions to these general features in the cases of a^, an, and Ogg are limited to small 

regions with a < 1 and /3 < 1 (top-right corners). For the kurtosis a^g , the size of the top-right exception 
region is significantly larger and, moreover, another even larger bottom-right corner appears. 

It is worth remarking that, while there still exists a prevalence of (cos 2 1 !)) < | (lifted-tennis-ball effect) 
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FIG. 1. Density plots, for uniform spheres (k = |), of the steady-state values of (a) the temperature ratio 9 = T r /Tt, 
(b) the dimensionless noise intensity 7 [see Eq. (12.1611 ]. (c) the translational kurtosis cl£q [see Eq. (12.27al) ]. (d) 
the rotational kurtosis o® [see Eq. (I2.27all [. (e) the scalar correlation cumulant ajy* [see Eq. (12. 27b 1) 1. and (f) the 
orientational correlation cumulant Uqq'* [see Eq. (12. 27b 1) 1. The contour lines correspond to (a) 8 = 0.1, 0.2,..., 0.9, 
(b) 7 = 0.1, 0.2,..., 1.5, (c) a$ = -0.01, 0.00,..., 0.09, (d) a$ = 0.00, 0.01,..., 0.05, (e) = 0.00, 0.01,..., 0.09, 

and (f) agg = -0.07, -0.06,..., 0.00. 


for the unforced granular gas in the HCSj 29 i 30 ' 36 ' 44 a significantly larger region with (cos 2 i9) > | (“cannon¬ 
ball” effect) appears close to a, /3 < 1 and a second region with the same behavior appears in the opposite 
perfectly smooth limit (a < 1, /? > — 1). The counterpart of this second region is completely absent in the 
case of a heated granular gas considered here. This is closely related to the fact that, while the quasismooth 
limit /3 —>• —1 is singular in the HCS j 29 i 36 it becomes regular in the heated case. Taking carefully the limit 
/3 —> — 1 in the analytical expressions for 9, 7 , and a '^, one gets 2 £ 

61 — 45a + 6 a 2 (l — a) 1 + /? n 2 \ 61 — 45a + 6 a 2 (l — a) 

15 — 11a + 2a 2 (l — a) 8 k + 1 ’ 7 1 ; 241 - 177a + 30a 2 (l - a) ’ 


(3.7a) 
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FIG. 2. Evolution, for a = 0.9, /3 = 0, and 7 ( 0 ) = — 6.381 of (a) temperature ratio 0(t) and reduced noise 

intensity (normalized to its initial value) 7 (t)/7 (0), and (b) cumulants o^(r). Time is measured in units of number 
of collisions per particle, r. Solid lines stand for DSMC results, dashed lines for the Sonine approximation, and 
dotted lines [only in panel (a)] for the Maxwellian approximation. 


-*20 


16(1-a) 


1 - 2 a 2 


241 — 177a + 30a 2 (l — a)'' 


x (0) 

a 02 


,( 0 ) 


z (1) 


0 . 


(3.7b) 


The expression of aJ,o' in Eq. (I3.7bl) coincides with the one derived directly for strict smooth spheres (/3 = 
— l)4i Interestingly, in the latter case the rotational distribution function (and hence the temperature T r ) 
is not uniquely defined since it preserves with time its initial arbitrary form. On the other hand, the 
limit P —> — 1 shows that the the rotational degrees of freedom become decoupled from the translational 
ones (a^' -7- 0, -7- 0) but the rotational distribution function tends to a Maxwellian -7- 0 with a 

temperature much smaller than the translational one ( T r /T t —> 0). 

Of course the theoretical results presented in this section are valid only if they agree with a direct compar¬ 
ison with a solution of the Boltzmann equation m free from any additional hypotheses. This comparison 
will be done in Sec. hyi 


IV. COMPARISON OF THEORETICAL RESULTS WITH DSMC DATA 

We display in this section detailed comparisons of the theoretical results with accurate DSMC data, as 
obtained from a code we wrote specifically for the Boltzmann equation m for this system. For data 
noise reduction, we have used 2 x 10 6 simulation particles for all cases in this work. More details on the 
implementation of DSMC in a heated granular gas can be found in Ref. [42j. We fixed the value k = | 
(uniform spheres) for the reduced moment of inertia and varied either the coefficient of normal restitution 
a (at fixed f3) or the coefficient of tangential restitution ft (at fixed a). The noise intensity parameter 
was typically chosen as Xo = 7 ywrcr 2 pt( 0 )/TO] 3//2 , which corresponds to 7 ( 0 ) = — 6.381, but other 

values were also considered to check the independence of the steady state on the precise value of Xo- hi all 
cases, the particle velocities were initially sampled from a Maxwellian distribution [i.e., 0) = 0] with a 

common temperature T r (0) = T t ( 0) [i.e., 0(0) = 1]. For each different system, the corresponding data for 
both transient and steady states were recorded. 


A. Transient behavior 

As an illustration of the relaxation towards the steady state, and for the case a = 0.9, /3 = 0, and 
7(0) = 77 — 6.381, we show in Fig. [2] the evolution of the temperature ratio 6(t) and the dimensionless 

noise intensity y(r) in panel (a), and of the four cumulants c 4 a^fr), a^(r), and Oqq \t) in panel 
(b). A very good agreement between theory and simulation is observed. In the case of 7 (r), the Maxwellian 
and Sonine approximations give results practically indistinguishable from the DSMC ones. For 0(r), the 
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(a) 


(b) 


(c) 



FIG. 3. Steady-state values of (a) the temperature ratio 9, (b) the reduced noise intensity 7 , and (c) the cumulants 
as functions of ft for a = 0.9. Circles stand for DSMC data [with 7 ( 0 ) = — 6.381], solid lines for the 

Sonine approximation, and dashed lines [only in panels (a) and (b)[ for the Maxwellian approximation. The DSMC 
data for the quantities b (x) and h (+) [see Eq. (12.281 1 are also included in panel (c). 


Sonine approximation is still virtually perfect, while the Maxwellian one tends to slightly overestimate this 
quantity. The initial reduced white noise intensity 7 ( 0 ) = is larger than the initial cooling rate a4o*(0), 

and therefore the external-force heating effect dominates over the inelastic cooling. As a consequence, 7 (r) 
monotonically decays in time [see Eq. (12.191) ] until the steady state is reached (after about 10 collisions per 
particle). In this respect, note that, according to Eq. (12.161) . the normalized quantity 7 (t)/ 7 ( 0 ) plotted 
in Fig. [2](a) is equivalent to [T i (r)/T t (0)] _3 / 2 . Thus, the steady-state value 7(t)/7(0) —>• 0.099 implies 
T t (r)/T t { 0) — > 4.66. The heating is much less efficient in the case of T r since the stochastic force acts 
directly on the translational velocity only. Therefore, the rotational-to-translational temperature ratio 
decays (but non-monotonically) to T r (r)/T t (r) 0.217 and the rotational temperature has hardly increased 
to T r {r)/T r ( 0) 1.010. 

The evolution of the cumulants is much less intuitive. As can be observed in Fig. [2jb) , for the case 
under consideration (a = 0.9, ft = 0) the translational kurtosis hardly increases, the rotational kurtosis 
a,Q 2 grows non-monotonically, and the more relevant non-Maxwellian properties are associated with the 
correlation quantities and cIqq'. All these features are well captured by our Sonine approximation, 
although it tends to slightly overestimate a^ and underestimate and |<* 00 ^ I * 


B. Steady state 

Now we focus on the steady state. The influence of the coefficient of tangential restitution /? on 0, 7 , and 
dj2 at a fixed coefficient of normal restitution a = 0.9 is shown in Fig. [3] The DSMC data for the correlation 
quantities b and h [see Eq. (12.281) ] are also plotted in panel (c). Panels (a) and (b) confirm the excellent 
performance of the simple Maxwellian approximation (13.21) and of the more elaborate Sonine approximation 
in what concerns 9 and 7 . The temperature ratio 9 monotonically increases from 9 —> 0 in the smooth limit 
/? —> — 1 [see (13.7aD ] to 9 ~ 1 at the perfectly rough case /? = 1. Although not visible on the scale of Fig. 
0]a) , the DSMC value at /? = 1 is actually 9 = 0.998. The quantity 7 exhibits a non-monotonic dependence 
on /3, as already discussed in connection with Fig. [Q[b), with a maximum 7 max — 0.675 at /3 ~ 0.3. 

Regarding the cumulants, the behaviors observed in Fig. ESA) agree with what might have been anticipated 
by imagining a vertical slice a = 0.9 in Figs, [ljc)— (f). While a^, a^, and |(Iqq | vanish at (3 —> —1, are 
very small at /3 = 1 , and present maxima in the medium-roughness region (3 ~ 0, the translational kurtosis 
has a more complex behavior with a (positive) maximum at /3 ~ 0.34 and two (negative) local minima 
at /3 ~ —0.79 and /? ~ 0.93. We can also observe that the expectations (13.61) are indeed satisfied. Actually, 
the only visible differences between a qq\ b, and h occur in the region /3 ~ 0, where one has \h\ < |agg*| < |6|. 
As is apparent from Fig. [31(c), the performance of the Sonine approximation is very good, although small 
discrepancies, especially in the case of the rotational kurtosis ajyj , appear in the region of medium roughness 
(/3 ~ 0), where the magnitudes of the cumulants are higher. Although hardly visible in Fig. (3](c), the theory 
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FIG. 4. Steady-state values of (a) the temperature ratio 9, (b) the reduced noise intensity 7 , and (c) the cumulants 
as functions of a for p — —0.9. Circles stand for DSMC data [with 7 ( 0 ) = — 6.381], solid lines for the 

Sonine approximation, and dashed lines [only in panels (a) and (b)] for the Maxwellian approximation. The DSMC 
data for the quantities b (x) and h (+) [see Eq. (12.281 1 are also included in panel (c). 




FIG. 5. Steady-state values of (a) the temperature ratio 9, (b) the reduced noise intensity 7 , and (c) the cumulants 
as functions of a for /? = 0. Circles stand for DSMC data [with 7 ( 0 ) = §\/ 7 r, f \/ 7 r, §\/ 7 r, 3y/n, ^y^r], solid 

lines for the Sonine approximation, and dashed lines [only in panels (a) and (b)] for the Maxwellian approximation. 
The DSMC data for the quantities b (x) and h (+) [see Eq. (12.281 ] are also included in panel (c). 


succeeds in predicting that a^.J < 0 (leptokurtic rotational velocity distribution), < 0 (high c correlated 
with low w, and viceversa), and clqq > 0 (cannon-ball effect) in the regions 0.89 < /3 < 1, 0.92 < /3 < 1, and 

0.94 < j3 < 1, respectively. Partially due to the fact that the translational kurtosis a:^ is the cumulant of 
smaller magnitude in the central region, its Sonine prediction is excellent for all (3. 

Now we turn to the analysis of the influence of a at fixed /?. To that end, we have chosen the characteristic 
values (3 = —0.9 (weak roughness), (3 = 0 (medium roughness), and j3 = 0.9 (strong roughness). The results 
are displayed in Figs. Q1 El and [HI respectively. Note that in Fig. [5] the plotted DSMC data correspond to six 
different values of the initial reduced noise intensity ranging from 7 ( 0 ) = — 0.532 to 7 ( 0 ) = Tpy^ — 

6.381. The almost perfect collapse of the points confirms the independence of the steady state on the values 
of Xo an d, equivalently, 7 ( 0 ). 

From Figs. Oja), [3(a), and[H|a) we observe that, even though the Maxwellian approximation states that 
the temperature ratio 6 is independent of a, this quantity tends to increase with increasing inelasticity, 
except near the elastic limit (a < 1) for /3 = —0.9 and 0.9, were a minimum value is reached. These fine¬ 
grained phenomena are very well accounted for by our Sonine approximation, although the latter tends to 
slightly overestimate 9 (note the magnified vertical scales). At fixed f3 a decrease of a produces an increase 
of dissipation and, consequently, an increase of 7 . According to Figs. Qlb), [3(b), and [3(b), the parabolic 
dependence of 7 on a predicted by the Maxwellian approximation [see Eq. (13.2[) ] tends to lie slightly below 
the DSMC points for small a. This deviation is satisfactorily corrected by the Sonine approximation. In 
what concerns the cumulants, Figs. Sic), [5](c), andilc) show again a very good predictive power of the 
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FIG. 6 . Steady-state values of (a) the temperature ratio 9, (b) the reduced noise intensity 7 , and (c) the cumulants 
a:^. as functions of a for (3 = 0.9. Circles stand for DSMC data [with 7 ( 0 ) = ^yTr — 6.381], solid lines for the 
Sonine approximation, and dashed lines [only in panels (a) and (b)[ for the Maxwellian approximation. The DSMC 
data for the quantities b (x) and h (+) [see Eq. (12.281) 1 are also included in panel (c). 


Sonine approximation. As in Fig. c), the larger discrepancies are found for for medium roughness 
(j3 = 0). Interestingly, the translational kurtosis becomes larger than the rotational one, i.e., > a^J 

if a is small enough. Also, the practical equivalence between the correlation quantities a ( yj. 6 , and h is 
maintained for all a, although small deviations can be observed at (3 = —0.9 and (3 = 0. 


C. Marginal velocity distribution functions 

The four fourth-degree cumulants , a^, and encapsulate the basic piece of information about 
the deviations of the true velocity distribution function </>(c, w) from the (two-temperature) Maxwellian 
4>m(c, w ) = 7r~ 3 exp(— c 2 — w 2 ). On the other hand, it seems interesting to perform a more direct comparison 
between <j)( c, w) and w ) and assess to what extent w) can be represented by the truncated Sonine 

expansion (IQ) . Given that (/)( c,w) depends on the three quantities c 2 , w 2 , and (c • w) 2 (or, analogously, 
cos 2 $), a plot of w), or of the ratio w)/((im{c, w ), would require a four-dimensional space. Therefore, 
it is convenient to analyze instead the projections of (f>( c, w) represented by the five marginal distributions 

H233D . 

Figure[7]displays the ratios R c (c) = (j> c (c)/<j> CtM (c), Rw(w) = (t> w (w)/<j> WlM (w), R cw (c 2 w 2 ) = </> C i I ;(c 2 w; 2 )/(/>c U ,,m(c 2 w 2 ), 
R c . w ((c • w) 2 ) = </> c . w ((c • w) 2 )/^ c . Wj m((c • w) 2 ), and Rg(cos 2 0) = 4>a(cos 2 $)/<A#,m(cos 2 0) for the rep¬ 
resentative system a = 0.9, (3 = 0. As we saw from Fig. 0(c) , at this value of the coefficient of normal 
restitution (3 the magnitudes of all the cumulants, except the translational kurtosis a^/ are close to their 
maxima for a = 0.9. As can be observed from Fig. [7](a), although the Sonine approximation predicts very 
accurately the small value of a.^j (a^ = 0.00144 versus the DSMC value 0.00137), the function R c (c) is only 
qualitatively captured by Eq. (I3.5a[) . In particular, the high-velocity tail of R c (c) (for c > 2.5) is higher than 
predicted by the Sonine approximation. In contrast, even though the Sonine approximation overestimates 
the rotational kurtosis by about 12% (ajyj = 0.0343 versus the DSMC value 0.0305), Fig. Qb) shows that 
the predicted function R w (w ) is rather accurate up to w ~ 3.4, where R w (w) reaches values as high as 2.5. 

Again, however, the high-velocity tail of R w (w) is more overpopulated than the Sonine prediction. In fact, 
the high-velocity tails of the velocity distribution, which typically do not influence the first few cumulants, 
are not expected to be accounted for by any truncated Sonine expansion. 

Now, we turn to the three marginal distribution functions related to the translational-rotational correla¬ 
tions. FigurejTJc) shows that small ( c 2 w 2 < 0.56) and large ( c 2 w 2 > 6.15) values of c 2 w 2 are more probable 
than it might be expected from the Maxwellian distribution. This phenomenon is well described by the 
Sonine approximation (I3.5bl) . although some differences are observed around the minimum of R cw (c 2 w 2 ). 

Anyway, the error in (c 2 w 2 ) — | is less than 3% (|o^ = 0.162 versus the DSMC value 0.167). According to 
Fig. EH), most of the values of (c • w) 2 within the “thermal” range are less probable than the Maxwellian 
expectations. This agrees with a negative value of ((c • w) 2 ) — |, which is very well described by the Sonine 
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FIG. 7. Ratios of marginal distribution functions (a) R c (c) = <t> c {c)/<t>c,M{c), (b) R w (w) = ij> w {w) / (/)w,m(w), (c) 
Rcw{c 2 w 2 ) = 0 C u;(c 2 w 2 )/(/)c«,,m(c 2 w 2 ), (d) i? c w((c ■ w) 2 ) = 0 c -w((c ■ w) 2 )/^ C w,m((c • w) 2 ), and (e) R&( cos 2 rl) = 
(f)$(cos 2 i?)/0^ : m(cos 2 i?) for a = 0.9 and /3 = 0. Symbols stand for DSMC results and lines for the Sonine approxi¬ 
mation theory [see Eqs. dSJ]. 


approximation (|a^ + = “0-0518 versus the DSMC value —0.0520). Note that the positive value of 

0 $ (i.e., ( c 2 w 2 ) > ( c 2 )(w 2 )) is dominated by a negative value of (i.e., 3((c-w) 2 ) < ( c 2 w 2 )), resulting in 
3((c-w) 2 ) < ( c 2 )(w 2 ) < ( c 2 w 2 ). Obviously, in order to preserve normalization, i? c . w (( c-w) 2 ) must be larger 
than unity for values of (c-w) 2 higher than those considered in our DSMC simulations (according to the So¬ 
nine approximation, this would happen for (c • w) 2 > 19.26), but they do not affect the sign of ((c • w) 2 ) — |. 
Finally, Fig.[7](e) confirms the lifted-tennis-ball effect, i.e., the values of cos 2 $ smaller than about 4 are more 
probable than one could expect from a Maxwellian distribution, so that (cos 2 t?) — ^ = < 0. The Sonine 

approximation predicts the former difference within an error smaller than 5% (j^b = —0.0169 versus the 
DSMC value —0.0177). Apart from the average value (cos 2 $), the Sonine approximation (13.5dl) describes 
rather accurately the whole distribution R#( cos 2 $), as Fig. [7][e) shows. 


V. CONCLUSION 

We have developed in this work a kinetic theory for a dilute granular gas of rough hard spheres that is 
heated by a uniform stochastic volume force. Our theory is based on a Sonine polynomial expansion of 
the granular gas distribution function around the Maxwellian. Thus, we obtain a numerical solution of the 
velocity distribution function at all times and an analytical solution for the uniform steady state that is 
reached after a relatively short relaxation time. We have tested our theoretical approach against DSMC 
data, obtaining a very good agreement between simulation and theory. 
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Contrary to the results for a freely cooling granular gas of rough spheres in a recent work*2S the cumulants 
of the distribution function reach low values (always below 0 . 1 , as measured from simulations and theory) 
and thus our theoretical approach works in the whole range of inelasticity and roughness values. In other 
words, it would not be necessary for the heated granular gas to add higher order terms in the truncated 
Sonine expansion that we have considered, as done in recent works on the HCS for inelastic rough spheres^ 
Also, thanks to the existence of a true steady state (where inelastic cooling and external heating balance each 
other), the quasismooth limit /3 —>• — 1 is regular and thus the results for purely smooth spheres (/3 = — 1 ) 
are recovered in that limit. Again, this differs from the situation in the HCSi 29 i 36 ’ 44 

Another relevant result for the heated granular gas is that the translations and rotations in particles, al¬ 
though still correlated as in the HCS^ show a very mince cannon-ball effect region and this effect completely 
disappears from the nearly-smooth region, appearing only near the completely elastic and rough limit. This 
is in contrast with the behavior in the HCS, where two neat cannon-ball regions may be found ; 29 ' 30 ! 44 

Summarizing, the white noise force has an effect of making the distribution function more similar to 
the Maxwellian for the whole range of coefficients of normal and tangential restitution. Since our Sonine 
expansion is done around the Maxwellian we expect the theoretical solution that we provide in this work 
to be an accurate description and a reference for researchers interested in the description of the transport 
properties of a heated granular gas of rough spheres. 
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Appendix: Collisional moments 


We write in this Appendix the Sonine-approximation expressions of the collisional moments (with 

p + q + 2r = 2 and 4) as linear functions of the cumulants (with j + k + 21 = 2) and nonlinear functions 
of the temperature ratio 6. Those expressions, where we we make use of the parameters (ESI) , are 
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